Rotational Stabilization of Magnetically Collimated Jets 
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ABSTRACT 

We investigate the launching and stability of extragalactic jets through nonlinear magnetohy- 
drodynamic (MHD) simulation and linear eigenmode analysis. In the simulations of jet evolution, 
a small-scale equilibrium magnetic arcade is twisted by a differentially rotating accretion disk. 
These simulations produce a collimated outflow which is unstable to the current driven m — 1 
kink mode for low rotational velocities of the accretion disk relative to the Alfven speed of the 
coronal plasma. The growth rate of the kink mode in the jet is shown to be inversely related to 
the rotation rate of the disk, and the jet is stable for high rotation rates. Linear MHD calcu- 
lations investigate the effect of rigid rotation on the kink mode in a cylindrical plasma. These 
calculations show that the Coriolis force distorts the m = 1 kink eigenmode and stabilizes it at 
rotation frequencies such that the rotation period is longer than a few Alfven times. 

Subject headings: galaxies: jets, galaxies: magnetic fields, plasmas, MHD, instabilities 
Introduction 



Large-scale, highly collimated energetic plasma 
outflows are observed in some active galactic nu- 
clei (AGN). Many models have been proposed for 



Many of the earlier computational efforts to 
model extragalactic jets concentrate on two- 
dimensional MHD models in which the accre- 



tion disk is treated as a boundary condition (Ro- 



the formation of these jets (Ferrari 1998), but 
their launching, collimation, and stability remain 
open issues. Recent observations indicate that the 
magnetic field structure in AGN jets is helical in 
nature ( |Asada |[2005| |Gabuzda ] 2004| |Marscher| 
2008 1 . This suggests that magnetic fields play 
a strong role in the collimation of AGN jets, as 
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Even though each of the studies cited uses a dif- 
ferent initial magnetic field, they all observe the 
formation of a steady outflow. More recent three- 
dimensional MHD simulations study the stability 



of the jet far from the galactic nucleus (Nakamura 



was proposed by Blandford and Payne (Blandford 



1982), and that one can use a magnetohydrody- 



namic (MHD) model to describe their formation 
and evolution. However, both theory and labora- 
tory experiments show that helical MHD equilib- 
ria can be unstable to current-driven kink modes. 
Understanding the effect of the kink mode on jet 
morphology is therefore critical to understanding 
their evolution. Here, we describe a computational 
MHD study of the stability of plasma jets relative 
to the kink mode and the effect that jet rotation 
has on the stability properties. 
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2001 ). These calculations inject flow and torsional 
Alfven waves into an MHD equilibrium and show 
that wiggled structures form in the jet due to the 
kink mode. Similar calculations, which consider a 
more realistic atmosphere into which the jet ex- 
pands, also examine the effect of the kink mode 
on the jet (Nakamura 2004). These calculations 



show that rapid rotation of the jet can have a sta- 
bilizing effect. The study discussed here aims to 
further examine the effect of equilibrium rotation 
on the stability of an expanding jet. 

The effect of equilibrium flow on current-driven 
MHD instabilities has been investigated both in 
theory and laboratory experiments. Linear MHD 
calculations show that a sheared axial flow has a 
stabilizing effect on the kink mode, while a uni- 
form axial flow has no effect on the growth of 
the instability (Shumlak 1995). This effect was 
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confirmed experimentally (Shumlak 2003). Later 



theoretical work studied the effect of sheared he- 
lical flow on the kink mode and showed that the 
sheared azimuthal flow stabilizes the mode by cre- 
ating a phase shift in the plasma eigenfunctions 
flWanex |[2005| [2007 1 . 

The work discussed here extends the two- 



dimensional simulations of jet launching ( 
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to three-dimensions via nonlinear MHD calcula- 
tions and considers the effect of jet rotation on the 
current-driven kink mode. By scanning the rota- 
tion of the disk, we scan jet rotation, and similar 
to previous results (Nakamura 2004), the rota- 



tion of the jet is observed to stabilize the column. 
To better understand the stabilizing mechanism 
of the rotation, we perform linear MHD analysis 
for a simple cylindrical plasma equilibrium with 
rigid rotation. These calculations show that the 
Coriolis force stabilizes the non-resonant kink by 
distorting the eigenmode. 

The paper is organized as follows. Section [2] dis- 
cusses the results of nonlinear simulations of ex- 
tragalactic jet launching and evolution. The sta- 
bility with regard to the kink mode is shown to 
depend on the rotational velocity of the accretion 
disk relative to the Alfven speed of the initial mag- 
netic arcade. Motivated by this result, Section [3] 
examines the linear stability of the kink mode in 
a cylindrical equilibrium with rigid rotation via 
initial- value MHD calculations. The results show 
that rigid rotation provides a stabilizing effect. In 
Section [3] ideal MHD eigenvalue calculations are 
used to confirm the results of Section [3] and to 
examine the effect of equilibrium rigid rotation on 
the unstable range of axial wave numbers. We also 
examine the physical mechanism of rotational sta- 
bilization using the eigenvalue calculations in Sec- 
tion 2] and show that the Coriolis force stabilizes 
the kink mode. Discussion of the results and con- 
clusions are given in Section [5j 

2. Nonlinear Calculations of Jet Propoga- 
tion 

To investigate jet propagation, we model the 
expansion of a magnetic arcade due to accretion 
disk rotation using a non-relativistic MHD model 
which ignores gravitational effects. Similar to pre- 



Ustyugova~||2000 1, the accretion disk is treated as 
a boundary condition on the computational do- 
main. The simulation is initialized with axisym- 
metric vacuum magnetic field that is tied to the 
disk and has zero net magnetic flux through the 
disk. Thus, both ends of all magnetic field lines 
are anchored to the accretion disk. The differential 
rotation of the accretion disk, which rotates with 
a Keplerian velicity profile, injects magnetic helic- 
ity and magnetic pressure into the magnetic field, 
causing it to coil and expand. The coiled mag- 
netic field produces a hoop stress on the plasma 
that collimates it on the central axis. The effect 
of jet rotation on the stability of the column is ex- 
plored by varying the rotation rate of the accretion 
disk in individual simulations. 

We numerically evolve the visco-resitive nonrel- 
ativistic MHD equations, 
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vious studies (Romanova 1997 Ouyed 1997 
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where n is the particle density, B is the magnetic 
field, v is the flow velocity, p is the thermal pres- 
sure, T is the ion and electron temperature, K is 
the isotropic thermal diffusivity, v is the viscosity, 
rj is the resistivity, and 7 is the ratio of the spe- 
cific heats chosen such that 7 = 5/3. The particle 
density, n, is related to the mass density, p, by a 
factor of the ion mass. The pressure and tempera- 
ture are related by the ideal gas relation, assuming 
that the electrons and ions have the same temper- 
ature, p = 2nksT. There is an extra term added 
to the right hand side of the continuity equation 
(Eq. [IJ, given by V-DVn. This diffusive term is 
added for numerical smoothing and the diffusivity 
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coefficient, D, is generally chosen to be small. The 
thermal diffusivity, K, is chosen to be 100 times 
the electromagnetic diffusivity. The effect of the 
gravitational force due to the massive galactic cen- 
tral object has been ignored, so gravity does not 
appear in the momentum equation (Eq. |3|. 

The MHD equations are evolved in time using 
the NIMROD code (ISovinec 112004b . NIMROD is 



well bcnchmarked and has been used to model a 
wide array of plasma experiments ( Sovinec ||2003 1 
and magnetospheric physics ( Zhu |2006 1 . A cylin- 
drical computational domain with a cylindrical co- 
ordinate system given by (r, 6, z) is used. The spa- 
tial discretization scheme combines two numerical 
methods. A mesh of high-order finite elements 
is used in the poloidal (r-z) plane, where the de- 
gree of the polynomial basis functions is chosen by 
the user, and the azimuthal (6) direction is repre- 
sented with finite Fourier scries. The parameter m 
is used to identify Fourier components in the az- 
imuthal direction. Convergence studies show that 
a resolution of < m < 5, is sufficient for the 
dynamics of the expanding jet. Using logarithmic 
packing of the poloidal mesh on the central axis 
and disk boundary, we can resolve the jet dynam- 
ics in a large domain using a poloidal mesh of 48 
by 48 fifth-order elements. 

Previous studies searching for steady state out- 
flows have treated the outer boundaries of the do- 
main with open boundary conditions, allowing ki- 
netic and magnetic energy to flow out of the do- 
main ( Romanova ||1997 Ouyed |[l997 Ustyugova 
2000 ). We use closed, perfectly conducting bound- 



ary conditions on the outer boundaries to avoid 
inward propagating wave characteristics. While 
this boundary condition is certainly unphysical, 
the outer boundaries are placed at a distance of 
r = z = 100 r», where r$ is the inner radius of 
the accretion disk, which is far from the dynamic 
region of the calculation. 

The model of the accretion disk/jet system 
treats the accretion disk as a boundary condition 
at 2 = 0, where a smoothed axisymmetric Keple- 
rian velocity profile is applied to v§: 



ve(r, 9,z = 0) 



\[GM r 
(r 2 + rf )3/4 



(5) 



The remaining components of the fluid velocity on 
the disk boundary at z — are chosen such that 



v r = v z = 0. A Dirichlet boundary condition is 
applied on the other, distant boundaries with v = 
0. On all of the boundaries, the number density 
is constrained to be constant. Mass is allowed to 
diffuse through the disk boundary (to fill in the 
coronal mass that is removed by the jet flow) by 
shaping the diffusivity parameter, D, in Eq. [TJsuch 
that it is large near the disk boundary and small 
in the rest of the domain. All of the boundaries 
are treated as perfect conductors by holding the 
normal component of the magnetic field constant 
in time. 

The initial condition is a currentless coronal 
magnetic field, just above the accretion disk. This 
field is chosen such that there is zero net mag- 
netic flux through the accretion disk boundary. 
The poloidal magnetic flux, -0, defined by B = 
VV'XVtf, is chosen to be 
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on the disk boundary, where a is a parameter with 
< a < 1. The initial magnetic field is found in 
terms of the magnetic potential, <&, which satisfies 
Laplace's equation and is given by 



B = V$. 



(7) 



An analytic solution for <I> in the domain is found 
by solving the boundary value problem V 2 $ = 0, 
where the normal component of V$ on the accre- 
tion disk boundary is specified by Eq. [6] 
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dz 



1 di/j 
r dr ' 
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and $ = on all of the other boundaries. The 
poloidal flux on the disk boundary increases from 
zero at r — to a maximum value at the O-point of 
the magnetic field, defined to be the radius where 
B z = 0, and then exponentially decays to zero. 
For the calculations discussed here, the values a = 
3/4 and r a = 30 Vi are used. For this choice of a, 
the radius of the O-point of the magnetic field is 
at 15.13 r,. 

The initial number density and temperature are 
constant across the computational domain. Thus, 
the entire domain is initially filled with a plasma 
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that is essentially unmagnetized away from the ini- 
tial arcade, and the magnetized jet expands into 
this thermal plasma. The initial flow velocity is 
set to zero, and the accretion disk flow is ramped 
from zero at t = to a steady profile within one 
turn of the disk at r — Vi. The Keplerian flow of 
the disk acts to twist the coronal magnetic field, 
building magnetic pressure above the disk which 
launches the outflow. This twisting of the mag- 
netic field also creates a strong 0-component to 
the field, causing a hoop stress which pinches the 
plasma on the central axis and collimates the out- 
flow. 

The results discussed here are given in units 
of the initial field quantities. All velocities are 
given in units of the Alfven velocity at the origin, 

— 1/2 

VA = B Q p , where B a and p Q are the magnetic 
field and mass density at the origin respectively. 
The magnetic field is given in units of B a . Time 
is given in units of Tj, the rotation period of the 
disk at r = T{. 

Four dimensionless parameters are used to de- 
scribe the system. The first three are commonly 
used to describe plasma systems: the Lundquist 
number, S = tr t^ 1 , where tr = fj, n rf rj^ 1 is 
the resistive diffusion time across the inner radius 
of the accretion disk and ta — r % "A^ 1 is the 
Alfvenic propagation time across the inner radius 
of the disk based on the Alfven speed at the origin; 
the magnetic Prandtl number, Pjyj = v fi a i]^ 1 ; 
and the plasma beta at the origin, = P T Pg , 
where Pt is the thermal pressure and Pb is the 
magnetic pressure. The last dimensionless param- 
eter, the drive parameter, Vd, is defined as 
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where va is the Alfven velocity. This parameter 
can be understood as the ratio of how fast the ac- 
cretion disk twists coronal magnetic field lines to 
how fast the information of this twisting propa- 
gates through the corona. In order to maintain 
a constant resistive diffusion time relative to the 
rotation period of the accretion disk in different 
simulations, the parameter S-Vd is held constant 
as Vd is varied. Three sets of parameters are con- 
sidered: Pjyj, 0, and the product S-Vd are fixed at 
1, 1, and 2007T respectively, and Vd is varied with 
the values 0.5, 1.0, and 4.0. The values of Vd 
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Fig. 1. — Cross section of the z-component of the 
fluid velocity, for Vd = 0.5 and 4.0, at times t — 
65.6 and 121.7 respectively. Velocity is shown 
in units of va - Note that the domain extends to 
100 Ti in radius and height. 



are chosen to be similar to previous studies (Moll 
1 |2008l |Nakamura ||2004] |Ouyed |[2003| which con- 
sider sub-Alfvcnic disk rotation, and to extend the 
disk rotation to the previously unstudied super- 
Alfvcnic regime. 

The z-component of the fluid velocity for the 
Vd — 0.5 and 4.0 calculations at t = 65.6 and 
121.7 Ti respectively is shown in Fig. [T] While 
a collimated outflow is produced for both values 
of Vd, non-axisymmetric structure forms in the 
column for the Vd = 0.5 case, due to the presence 
of an MHD instability. The effect of the instability 
on the magnetic structure of the jet can be seen in 
Fig. [2] Here the magnitude of the magnetic field 
is shown for V D = 0.5 and 1.0 at t = 43 T l . While 
the jet has expanded to a similar length for both 
values of Vd at these times, the modification of the 
magnetic structure is more significant for Vd = 
0.5. For both cases an m = 1 kink mode creates a 
helical distortion to the magnetic structure. 

To confirm the source of the asymmetry in the 
Vd = 0.5 simulation, we plot the energy of indi- 
vidual Fourier components in Fig. [3j The m = 1 
component is the first to become unstable, and it 
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Vb = 0.5 V D = 1.0 




Fig. 2. — Three-dimensional contours of the mag- 
nitude of the magnetic field (B = 0.34 B ) for 
Vd — 0.5 and 1.0 after the m = 1 kink mode has 
saturated at t — 43 TJ. Note that the domain 
extends to 100 r-j in radius and height. 

nonlinearly drives the m > 1 components when it 
reaches a significant level. The nonlinear drive is 
confirmed by artificially resetting the dependent 
fields in the m = 1 component to zero during 
the course of a simulation. As can be seen from 
the dashed traces in Fig [3j removing the m = 1 
component causes the larger-m components to de- 
cay, until the m = 1 returns to a significant level. 
Thus, the m = 1 component nonlinearly channels 
energy into the m > 1 components. 

A plot of the magnetic energy of the m = 1 
Fourier component for all three jet simulations is 
shown in Fig. |4j The jet is unstable to an m = 
1 mode for Vd — 0.5 and 1.0, while it remains 
nearly stable for Vd — 4.0. We calculate the linear 
growth rate of the m = 1 mode by making a linear 
fit to the magnetic energy of the m = 1 mode when 
it is in the linearly growing phase. The growth 
rates for Vd = 0.5, 1.0, 4.0 are found to be 2.98, 
1.88, and 0.13 T" 1 respectively. Thus, we see that 
as the accretion disk rotation increases relative to 
the Alfven velocity of the coronal plasma, i.e. as 
the drive for the jet increases, the growth rate of 
the m = 1 kink mode decreases. 
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Fig. 3. — Solid lines show the energy of individ- 
ual azimuthal Fourier components as a function 
of time for Vd — 0.5. The dashed lines show re- 
sults from a second computation where the fields 
in the m = 1 component are reset to zero at 
t = 8.0 Tj. The magnetic energy is normalized to 
Vdo m Bo/(2fi ), where Vd om is the volume of the 
computational domain. 

3. Linear Initial Value Calculations 

The stability of the kink mode for large values 
of Vd indicates that flow plays an important role 
in stabilizing the jet. Motivated by this observa- 
tion, we perform linear initial value calculations in 
a simple geometry where the flow can be scanned 
systematically and consider both rotation and ax- 
ial flow. A cylindrical domain is used with a co- 
ordinate system given by (r, 6, z). The fields are 
defined to be periodic in the z-direction, and the 
boundary at r = r a is treated as a perfect conduc- 
tor, where r a is much smaller than the radius of 
the domain of the nonlinear simulations described 
in Sec. [2] The results are given in terms of the 
Alfven propagation time across the radius of the 
cy finder, ta = r a v^ 1 , where va is the Alfven 
speed at r = 0. Here, we solve a linear version 
of Eqns. [l||4]for perturbations to MHD equilibria, 
with an arbitrary perturbation included in the ini- 
tial velocity field. The dissipation coefficients are 
chosen to give a Lundquist number of S = 1 X 10 6 
and a magnetic Prandlt number of Pm = 1. If 
an MHD equilibrium is unstable, the solution ob- 
tained will be the most unstable linear eigenmode, 
and the growth rate is determined from the result- 
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Fig. 4. — Energy of the m = 1 Fourier mode as 
a function of time for Vd = 0.5 (solid line), 1.0 
(dashed line), and 4.0 (dashed-dotted line). For 
Vd = 0.5 and 1.0 the jet is observed to be unsta- 
ble to the m = 1 kink mode, while for Vd =4.0 
the jet is nearly stable. The magnetic energy is 
normalized to Vdo m B 2 /(2fj, ), where Vd om is the 
volume of the computational domain. 



ing exponential growth. 

Our MHD equilibria are based on the param- 
agnetic pinch ( |Bickerton ||1958| , which is a one- 
dimensional Ohmic equilibrium with uniform ax- 
ial electric field. The equilibrium is characterized 
by the parallel current profile, A(r), defined as 



^ Jo(r)-B (r) E 0z B 0z (r) 
A ( r ) = ET732 = „ D /^g » ( 10 ) 



^o(r) 2 



where a subscript is used to represent equilib- 
rium fields. The profile discussed here is defined 
in terms of the on-axis parallel current, A D = \{r = 
0), and the width of the equilibrium current pro- 
file decreases with increasing A Q . Given that the 
— 1/2 J A <5E* • SB (fx term is the only poten- 
tially destabilizing term in the linear ideal poten- 



tial energy that is independent of Vp Q ( Freidberg 



1987a), the parallel current is related to the free 



magnetic energy available to drive the kink mode. 
Moreover, for the paramagnetic pinch, A serves 
as a stability parameter for the mode. The equi- 
librium magnetic field is found by choosing a value 
for A and numerically integrating Ampere's Law, 
VxB = ^oJoi using Eq. [l0]for the parallel com- 
ponent of Jo- 



Fig. 5.— Radial A profiles for z = 20.25, 30.14, 
40.41, and 50.23 r,, from the nonlinear jet calcu- 
lation with Vd = 4.0, at t = 121.7 Tj are shown as 
solid colored lines. The curves overlap since there 
is not a significant change in A in the z-direction. 
The A profile for the paramagnetic pinch with 
A„ = 5.0 is shown as a dashed line. 



A plot of radial profiles of A from the nonlinear 
jet calculation with Vd — 4.0 at t = 121.7 Tj is 
shown in Fig. [5] for z = 20.25, 30.14, 40.41, and 
50.23 The curves overlap since there is not a 
significant gradient in A in the z-direction. Thus, 
a one dimensional equilibrium for the linear calcu- 
lations is a good approximation of the A profiles in 
the nonlinear jet calculations. For comparison, the 
A profile for the paramagnetic pinch with A Q = 5.0 
is also plotted in Fig. [5] 

The stability of diffuse pinches, such as the 
paramagnetic pinch, without equilibrium flow rel- 
ative to the ideal kink mode has been well stud- 
ied and is known to depend on the pitch of the 
magnetic field, P(r) = r B z (r) Bgir)" 1 . Con- 



sidering eigenfunctions of the form e 
ergy analysis shows that for m = 1 and 2?? = 0, 
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the plasma is stable if kP(r) > 1 or kP{r) < 
(k 2 r 2 - 1) (3 + fcV)" 1 for r > and all val- 
ues of k ( Robinson |1971 ) . When there is a region 
in the plasma where k'P(r) < 1 and k'P(r) > 
(k' 2 r 2 — 1) (3 + fc /2 r 2 ) -1 , there is a source of free 
energy for the m = 1, k = k' kink mode, and it 
may be unstable. When k'P(r) < 1 in the en- 
tire plasma, the mode is non-resonant. If there 
is a radius, r s , in the plasma where k'P(r s ) = 1, 
r s divides the plasma into two regions; one where 
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Magnetic pitch, P(r) = 
1 , for the paramagnetic pinch 



Fig. 6.— 

r B z (r) B e {r) 
equilibrium with A D = 5.0 and = 1.0 



there is free energy for the kink, and one where 
there is not; and the mode is called resonant. For 
the paramagnetic pinch, P(r) decreases monoton- 
ically, and there is free energy for the kink in the 
region with r > r s . Since the free energy for the 
kink is at radii larger than r s , the stabilizing effect 
of the conducting boundary at r — r a affects both 
resonant and non-resonant modes. The magnetic 
pitch profile of the equilibrium used here is shown 
in Fig. [6] Since P{r = 0) = 2.0 A~\ k > 0.5 A G 
modes are resonant and k < 0.5 A G modes are non- 
resonant. 

To examine the effect of jet rotation, we con- 
sider MHD equilibria with rigid rotation in the 
0-direction, and use £1 to denote the rotation fre- 
quency. While previous studies have shown that 
sheared flow is more efficient at stabilizing the 
kink mode (Wanex 2005), our nonlinear compu- 



tations show little azimuthal shear in the vicinity 
of the the jet. Radial profiles of the jet rotation 
frequency from the nonlinear jet simulation with 
Vd = 4.0 at various times and axial positions are 
shown in Fig. [7] As time increases, the jet ro- 
tation frequency reaches a steady state at higher 
axial positions along the length of the column. For 
all values of z, the rotation frequency is uniform 
to within 20% across the radius of the jet, which 
has a width of r < 1, and as the column propa- 
gates, the rotation frequency flattens. Thus, rigid 
rotation is a reasonable simplification. 

The paramagnetic pinch is often considered to 
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Fig. 7. — Radial profiles of the jet rotation fre- 
quency for z — 20.25, 33.44, and 50.23 r,, from the 
nonlinear jet calculation with Vd — 4.0, at times 
t = 80.6, 101.2, 121.7 Ti. The rotation frequency, 
O, is given by Q, = ve r~ l . 

be a force-free equilibrium in which the current 
is purely parallel to the magnetic field. However, 
equilibrium azimuthal flow breaks the force-free 
nature, since the centrifugal force of the flow must 
be balanced by another MHD force. Two choices 
of force balance are considered here. The first, 
labeled 'magnetic-balance', balances the centrifu- 
gal force against the force from the perpendicular 
current, 



J XB = -po^ 2 r, 



(11) 



and the parallel current is unchanged. Thus, while 
the current profile is modified by the introduction 
of the rotation, the A(r) profile, which is related 
to the free-energy source for the kink mode, is un- 
affected. The second force-balance model, labeled 
'pressure-balance', balances the centrifugal force 
against the equilibrium pressure, 



Vpo = p Q n r. 



(12) 



For this case the current profile is unchanged by 
the introduction of the rotation. However, as the 
rotation increases, the pressure profile becomes in- 
creasingly hollow in the sense that it peaks on the 
edge of the plasma, which can have a stabilizing 
effect (Freidberg 1987a). The equilibrium pres- 



sure is characterized by the plasma (3 on the cen- 
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tral axis. The choice of the values for A„ and is 
motivated by our nonlinear jet calculations, giving 
A = 5.0 and = 1.0. 

Our numerically computed growth rate of the 
m = 1, k = 0.4 A kink mode as a function of 
equilibrium rotation frequency, for both force bal- 
ance models, is plotted in Fig. [8j The results 
show that the growth rate of the mode decreases 
as rotation increases for both force balance mod- 
els. The growth rate decreases somewhat faster in 
the pressure-balance model than in the magnetic- 
balance model, which we surmise is a result of the 
additional stabilizing effect of the hollow pressure 
profile in the pressure-balance model. While the 
results point to rotation as the important stabi- 
lizing mechanism, force-balance requires changes 
to the pressure profile or the perpendicular cur- 
rent profile as rotation is increased. To examine 
the effect of modifying the equilibrium forces to 
balance the centrifugal force from the rotation, a 
plasma which has the same equilibrium current 
as the magnetic-balance model, but without rota- 
tion, is considered. Here, the equilibrium pressure 
gradient replaces the centrifugal force by defining 
a profile which is peaked on the central axis. The 
resulting growth rate is also plotted in Fig. [8] As 
the pressure gradient increases, the growth rate of 
the kink mode increases. This result confirms that 
rotation is the stabilizing influence in the fi-scans. 

Previous theoretical and experimental studies 
show that sheared axial flow can stabilize the kink 



mode in a cylindrical plasma (Shumlak 1995 



2003). Thus, we consider what effect axial flow has 



on jet stability in the nonlinear simulations via lin- 
ear initial value calculations with equilibrium axial 
flow. Non-rotating force-free paramagnetic pinch 
equilibria with Gaussian axial flow profiles, given 
by v z {r) = vm er( 2r / w ^ , are considered. Moti- 
vated by the axial flow profiles in the nonlinear jet 
simulations, we choose vm = 0.3 va and consider 
a range of w g from 5.0 to 50.0 A" 1 , where smaller 
values of w g correspond to larger flow shear. Axial 
flow profiles from the stable Vjj — 4.0 jet simula- 
tion and the Gaussian profile used for the linear 
calculation with w g — 25.0 1 are shown in Fig. 
[9] Growth rates of the kink mode as a function 



of w g are plotted in Fig. 10 A flow shear range 



comparable to that considered by Shumlak ( 1995 1 
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Fig. 8. — Growth rates from linear initial 
value calculations for the magnetic-balance and 
pressure-balance models, and with an equilib- 
rium pressure profile which replaces the centrifugal 
force (peaked-pressure model). For the peaked- 
pressure model, there is no equilibrium rotation; 
instead, the rotationally equivalent pressure is 
given by po(^, r) = — fi 2 r 2 /2. 



ference in the equilibria considered here and that 
examined by Shumlak (1995 1. Based on these 



is considered, but the change in the kink growth 
rate is less than 6.0%. We attribute this to the dif- 



results, we conclude that axial flow does not sig- 
nificantly influence the stability of the magnetic 
column in our nonlinear jet simulations. 

4. Linear Eigenvalue Calculations 

The results of the linear initial value calcula- 
tions indicate that the nonlinear simulations re- 
main robust to the kink mode for high rotation 
rates of the accretion disk because of the rotation 
of the jet itself. To further examine the effect of 
azimuthal rotation on the kink mode, we investi- 
gate the linear ideal MHD spectrum for rotating 
paramagnetic equilibria. This cigenmode analy- 
sis helps us develop physical insight into the effect 
of rotation, which is difficult to obtain from the 
initial value calculations. The theory considers 
a cylindrical domain which is periodic in the z- 
direction with the one-dimensional rigid-rotation 
equilibria described in Section [3] 

4.1. Linear Eigenvalue Theory 

The simplest approach in considering an MHD 
equilibrium with flow is to work in a Lagrangian 
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0.35 





Fig. 9. — Axial flow profiles from the stable non- 
linear jet simulation with Vjj = 4.0 at t = 121.7 Tj 
at the indicated axial locations. The black dashed 
line shows the Gaussian flow profile used in the 
linear calculations with vm — 0.3 va and w g = 
25.0 Aq 1 . 

representation. Assuming perturbations to the 
equilibrium depend on time as e~' luJt , the linearized 
MHD equation of motion is given by 

-p uJ 2 € - 2ip a ujv -'V$, 
+PoVo-V(vo-VC)=F(C), (13) 

where from this point, fields without subscripts 
represent perturbations and are assumed to be 



small (Freiman 1960 Waelbroeck 19961. The 



Fig. 10. — Growth rates from linear initial value 
calculations with equilibrium axial flow as a func- 
tion of the Gaussian v z profile width. The black 
dashed line shows the growth rate without equi- 
librium flow. 

As a consistency check, we also evaluate the 
spectra derived from an Eulerian frame of refer- 
ence by including equilibrium flow in the definition 
of the Lagrangian displacement vector, £, which 
satisfies ( Chandrasekhar ||1961 1 



dt 



Vx(£xv ). 



(18) 



plasma displacement, £, is defined by 



Linearizing the MHD equations in the Eulerian 
frame with rigid equilibrium rotation gives the fol- 
lowing momentum equation, force operator, in- 
duction equation, and pressure equation respec- 
tively, 



The linear force operator, F(£), is given by 



(14) 



uj 2 D i - 2iQu D (zx$) + rfi 2 (V-£) 



mf2 \ 

3+ f 

u D ) 



Po 



(19) 



F(£) = -Vp+— J xB 



Ho 



H (VxB)xB 

Mo 

-V-(po^Vo-Vvo), (15) 



where the perturbed magnetic field and the per- 
turbed pressure are given by 

B = Vx(£xB ), (16) 
f> = -(€-Vpo + 7PoV-€). (17) 



F(€)= ^(B .VB + B.VB ) 



V (p 



-BBi 



B = Vx(£xB )- 



U1 D 



(20) 



(V-£)(rft0-mz), (21) 



P=-(^^-^[l + —j{V^), (22) 
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mfl is the Doppler shifted eigen- 



where u>d — 
frequency. 

Generalizing the analysis in Ref. 



( Freidberg 



1987b), the linearized equations in either refer- 



ence frame are reduced to a pair of coupled first- 
order differential equations for the radial plasma 
displacement, £ r , and the total perturbed plasma 
pressure, P = p + B B /p o . Assuming spatial 
dependence of the perturbed fields of the form 



e „.,„„ „,„^ Eqg [i3|i7| ari d Eqs. [T&||22l become sys- 
tems of ordinary differential equations (ODE's) 
with respect to the r-coordinate. By considering 



the projection of Eqs 
rections, where b 



13 



si 

I Bo I 



and 



19 



in the b and fj di- 
and r) = bxf, the b and 
fj components of the plasma displacement can be 
solved analytically. Substituting these results into 
Eqs. [T"3| and [T7| and Eqs. fl9] and [22] produces sets 
of coupled ODE's with the same general form, 



A(r,u>) 



dr 



P 



= B(r,u) 



P 



(23) 



in both reference frames. 

We consider the plasma to be surrounded by 
a conducting shell at the radius r = r a by defin- 
ing £, r ( r a) — 0- The regularity condition at r — 
is imposed by the cylindrical geometry of the do- 
main. Expansion of £ r in a power scries for small 
values of r shows that regular solutions satisfy 
£ r oc r 171-1 . Equation 23 coupled with these 
boundary conditions defines an eigenvalue prob- 
lem with lj as the eigenvalue. 

It should be noted that while the form of this 
eigenvalue equation is the same in both reference 
frames, the ODE coefficient matrices A and B are 
unique to each frame. Equation [23] is derived for 
a general equilibrium flow in a Lagrangian frame, 
and the coefficients can be found in Bondcson 
(1987). The ODE coefficients for a plasma equi- 



librium with rigid rotation and uniform axial flow 



in an Eulerian frame can be found in Appl ( 1992 1 



Due to the complexity of the ODE coefficients 
in Eq. [23j we use a shooting method to solve the 
eigenvalue problem. A value is chosen for uj, and 
Eq. 23 is numerically integrated from r = to 
r = r a using fourth-order Runge-Kutta integra- 
tion. The choice of lu is varied until the eigen- 
function satisfies fr(?"a) = 0. A Newton-Raphson 
method is used to search the w-parameter space 



Magnetic Force Balance 




Fig. 11. — Growth rates of the non-resonant m = 
1, k = 0.4 A Q kink mode as a function of the equi- 
librium rotation frequency, f2, from Lagrangian 
and Eulerian eigenvalue calculations, and from the 
linear initial value calculations. 

for functions that satisfy this boundary condition. 

In the absence of equilibrium flow, the MHD 
force operators in Eqs. |T"5| and [20| are self-adjoint, 
and ui is either purely real or purely imaginary 
( Freidberg |1987c" I . With the introduction of equi- 
librium flow, the force operator is no longer self- 



adjoint, and u> and £(r) can be complex (Freiman 



1960 ). The real component of the eigenvalue, 5f[w], 
gives the oscillation frequency of the eigenmode, 
and the imaginary component, determines 
its growth or decay rate. The Newton-Raphson 
method employed here is generalized to search the 
complex parameter space (Press 2007). While 



Newton-Raphson readily generalizes to multiple 
dimensions, it converges only if the initial guess 
for the root is in the vicinity of the actual root. 
Since lj is either purely real or purely imaginary 
without the equilibrium flow, Newton-Raphson is 
used in a one-dimensional space to find Q[u>] with 
= for a given mode. The O = result is 
then used as an initial guess for a nearby equilib- 
rium with flow, and that result is used as an initial 
guess for a slightly larger value of il. This process 
is repeated for increasing values of fl. 

4.2. Linear Eigenvalue Results 

Results of the eigenmode analysis and growth 
rates from the initial value formulation of Sec. [3] 
for the non-resonant m = 1, k = 0.4 A D kink mode 
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Fig. 12. — Growth rates of the resonant m = 1, 
k = 0.6 A kink mode as a function of the equilib- 
rium rotation frequency, £1. 

can be seen in Fig. [IT] Here, calculations are 
shown for both both force-balance models in both 
reference frames. The curves from the two refer- 
ence frames are indistinguishable in this plot, and 
comparison of the eigenvalue formulation and the 
initial value formulation of the problem are shown 
to be in agreement. These results show that as 
increases, the growth rate of the kink mode de- 
creases and is stable with sufficient rotation. We 
note that the marginal rotation period is larger 
than the Alfven propagation time, i.e. Alfvenic 
flow within the cylinder is not required for stabi- 
lization. 

We also examine the effect of rotation on res- 
onant kink modes via the eigenvalue formulation. 
Growth rates for the m = 1, k = 0.6 A D mode 
can be seen in Fig. |12| While equilibrium rota- 
tion fully stabilizes the non-resonant kink mode 
described previously, rotation only reduces the 
growth rate of the resonant mode and does not 
completely stabilize it. 

The eigenmode solutions treat the radial bound- 
ary at r = r a as a solid wall by setting £(r a ) = 0. 
However, there is no close boundary surround- 
ing the plasma column in the nonlinear jet sim- 
ulations. To evaluate the influence of the wall 
location, we recompute the eigenvalues as r a is 
varied. The critical rotation frequency, fi c , for 
stabilization of the m = 1, k = 0.4 A Q kink mode 
as a function of r a is plotted in Fig. [13] The re- 
sulting critical rotation frequency asymptotically 



Fig. 13. — Critical rotation frequency for stabi- 
lization of the kink mode (Solid Line) , and growth 
rate of the kink mode for = (Dotted Line) , as 
a function of the outer radial boundary, r a , for the 
m = 1, k = 0.4 A Q kink mode. 

approaches the value il c = 0.24 (k va)~ x , indi- 
cating that the stabilizing effect of the rotation 
remains as r a — > oo. The growth rate with 17 = 0, 
lo{fa), as a function of r a is also plotted in Fig. 
13 The f2 c (r a ) and 7 D (^ Q ) curves follow the same 
asymptotic trend, implying that the dependence 
of Q, c on r a is related to the free energy of the kink 
mode and not due to any changes in the stabilizing 
influence of rotation. 

The linear eigenvalue formulation allows for the 
examination of a range of axial wave numbers. 
The growth rate of the m = 1 kink mode as a func- 
tion of k for various values of f2 is calculated, and 
the results are shown in Fig. [14] Without equilib- 
rium rotation, there are lower and upper bounds 
on the unstable values of k. Both the upper and 
the lower bound increase with increasing rota- 
tion. For the equilibria considered here, modes 
with k < 0.5 A are non-resonant, and modes with 
k > 0.5 A Q are resonant. While the range in k- 
space of unstable non-resonant kink modes de- 
creases with increasing rotation, the range of un- 
stable resonant modes broadens with small growth 
rates on the order of 10~ 3 . 

We have explored a range of j3 values to exam- 
ine the effect of equilibrium thermal pressure on 
rotational stabilization. For moderate values of /3, 
rotational stabilization is observed to be indepen- 
dent of (3. However, for low /3-values ((3 < 0.06) 
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Fig. 14. — Growth rates of the m — 1 kink mode as 
a function of the axial wave number, k, for various 
equilibrium rotation frequencies. 

rigid rotation destabilizes the kink mode for f2 > 
0.4 t^ 1 . The destabilized modes are compressible 
with a ^-component of £ that is much larger than 
the other components. Thus, these modes are sta- 
bilized by equilibrium pressure for the moderate 
values of f3 relevant to extragalactic jet systems. 

The choice of initial and disk boundary con- 
ditions in simulations of jet formation can have a 
profound effect on the magnetic pitch profile, P(r) 
flMoll ||2008[ ). Thus far, we have considered only 
equilibria with monotonically decreasing P{r) as is 
observed in the nonlinear jet simulations discussed 
in Sec. [2] To check the effect of rotation on a 
monotonically increasing pitch profile we consider 
equilibria with P(r) = 1/2 + r 2 /2. The growth 
rate of the m = 1, k — 1.0 r^ 1 mode as a function 
of equilibrium rotation frequency is plotted in Fig. 
[15] for both force balance models. Similar to the 
decreasing P(r) cases, rigid rotation is shown to 
stabilize the kink mode, and we conclude that the 
rotational stabilization mechanism is not sensitive 
to the shape of the P(r) profile. 

We also use the cigenmode calculations to inves- 
tigate the physical mechanism for the rotational 
stabilization. The linearized momentum equation 
in the Eulerian frame is given by 



— + pov-Vvo + poVo-Vv + pvo-Vvo = F(v), 

(24) 



Fig. 15. — Growth rate of the m = 1, k = 
1.0 r^ 1 non-resonant kink mode for monotonically 
increasing magnetic pitch equilibria as a function 
of equilibrium rotation. 

and the growth rate of the m — l,k = 0.4 A Q 
kink mode is calculated as a function of fi, remov- 
ing one equilibrium flow term from the left side at 
a time. The results are plotted in Fig. [16] but 
results with the pvo-Vvo term removed are not 
shown, as this term does not have a significant ef- 
fect. In the computations without the pq'v-'Vvq 
term, the growth rate increases with increasing f2, 
so this term must play a central role in the sta- 
bilization. With rigid rotation, this inertial term 
is 

(v-V)vo = -i ojd(zx£) + ft 2 (V-£)r. (25) 

By individually removing each of the two terms on 
the right side of Eq. [25] at a time, we have deter- 
mined that it is the first term which provides the 
stabilization. This term contributes to the Coriolis 
force in the frame of the plasma. 

Plots of the £ r and components of the eigen- 
function for various equilibrium rotation rates are 
shown in Figs. [T"7]and[T8] While there is a slight 
change in 5R[£ r ] and as O is varied, the change 
in Q[£ r ] and 5R[£e] is more apparent. We note that 
the Coriolis term locally couples the radial and az- 
imuthal components of translation due to the kink. 
This distorts the mode giving a radially dependent 
phase shift in £ r , and a corresponding change in 
the real part of £g. This result is similar to that 
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0.5 



Fig. 16. — Growth rates of the m = 1, k = 0.4 A Q 
kink mode as a function of rotation frequency with 
individual inertial terms removed from the lin- 
earized momentum equation. For the solid curve 
all of the terms are present, for the dashed curve 
the /9oVi-Vv is removed, and for the dot-dashed 
curve the p v -Vvi is removed. 



described in Ref. ( Wanex |2005 1, where a radially 
dependent phase shift in the cigcnmodc due to a 
sheared equilibrium flow is shown to stabilize the 
kink mode. Here, we find that a rotational flow 
without shear introduces a stabilizing distortion 
of the mode via the Coriolis force. 

It should be noted that the Coriolis term also 
appears in the po v o*Vv term in Eq. [24] 

(vo-V)v = — i O ujd(zx£) +m f2 £ 

+(l 2 (V-£)r -imr ft 2 (V-£) 0, (26) 

but when the (vo*V)v term is removed, the sta- 
bilization effect is not lost. Equation [26] contains 
another term which is first order in f2 given by, 
m O lod This term provides the Dopplcr shift 
in the frequency uj. This Dopplcr shift appears 
in the other MHD equations as well. Thus, re- 
moving the jOo v o*Vv term temporally decouples 
the velocity field from the magnetic field, reduc- 
ing the growth rate of the instability, as shown in 
Fig. [Te] 

We also examine the effect of rotation on the 
resonant eigenmodes. Plots of £ r for the m = 1, 
k = 0.6 A Q kink mode, for various equilibrium ro- 
tation rates, are shown in Fig. |19| Similar to 





Fig. 17. — The £ r component of eigenfunctions 
of non-resonant m = 1, k = 0.4 A D kink modes 
for various equilibrium rotation rates. The eigen- 
modes are normalized to the maximum value of 



the non-resonant case, the rotation introduces a 
significant phase shift in the radial component of 
the cigcnfunction. However, for the resonant case, 
there is also a significant change in the real part 
of £ r near the rational surface. 

To assess the rotation in the simulated mag- 
netic columns described in Sec. [2] we calculate the 
rotation frequencies at different values of z. The 
rotation frequencies plotted in Fig. [20] are deter- 
mined by making linear fits to the ^-component of 
the fluid velocity over the radial coordinate. The 
simulation times chosen for these profiles are such 
that the kink mode is in the linear phase for the 
Vd = 0.5 and 1.0 calculations, as can be seen in 
Fig. [4] It is clear that angular momentum in- 
jected by the accretion disk is transported axially 
by the jet as it expands. As Vd increases, the ro- 
tation rate of the jet increases, providing greater 
stability for the kink mode. 

For comparison to the results of the linear MHD 
calculations, we examine the m = 1 kink mode in 
the nonlinear jet simulations when it is in the lin- 
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Fig. 18. — The £g component of eigcnfunctions 
of non-resonant m = 1, k = 0.4 A c kink modes 
for various equilibrium rotation rates. The eigen- 
modes are normalized to the maximum value of 

ear phase. The_m^= 1 Fourier component of v r is 
plotted in Fig. 



21 



for the unstable Vd = 0.5 and 
1.0 jet simulations at times t = 8.76 and 10.26 Tj, 
respectively. For these times, the kink mode is in 
its linearly growing phase. Since the modes plot- 
ted in Fig. [2l]extend across the entire width of the 
jet, we conclude that the kink mode observed in 
the nonlinear simulations is a non-resonant mode. 
According to our linear results, these modes would 
be stable with increased rotation, as is the case in 
the Vd — 4.0 simulation. Similar to the eigen- 
modes from the linear analysis shown in Fig. |17| 
the distortion of the linear eigenmodes in the jet 
simulations (Fig. 21) is due to a radially depen- 



dent phase shift in v r . 

5. Discussion and Conclusions 

Nonlinear non-relativistic MHD simulations of 
jet evolution, starting from an equilibrium coronal 
plasma with zero net magnetic flux through the 
accretion disk, show the formation of a collimated 
outflow. This outflow is unstable to the current 



Fig. 19. — The £ r component of eigcnfunctions of 
resonant m = 1, k = 0.6 A Q kink modes for various 
equilibrium rotation rates. The eigcnfunctions are 
normalized to the maximum value of 9?[£ r ]. The 
vertical black line shows the position of the ratio- 
nal surface at r = 0.232 r„. 



driven m = 1 kink mode for low rotation veloc- 
ities of the accretion disk relative to the Alfven 
speed of the coronal plasma. As it saturates, the 
kink mode broadens the outflow, but does not de- 
stroy the collimation. Similar to previous results 



(Nakamura 2004), for large rotation velocities of 
the accretion disk, the outflow is shown to be sta- 
ble against the kink mode. Moreover, the growth 
rate of the m = 1 kink mode is shown to be in- 
versely related to the rotation rate of the accretion 
disk. This result is counter-intuitive in the sense 
that as the accretion disk rotates faster, the colli- 
mating magnetic field in the jet coils tighter. As 
the coiling of the magnetic field increases, the cur- 
rent increases. Since the current is the source of 
free energy for the kink mode, one would expect 
that the jet would be more unstable for high ro- 
tation rates of the accretion disk. However, we 
observe that it is stable in this regime. 

Motivated by the result of the nonlinear jet sim- 
ulations, we explore the effect of rigid rotation on 
the m = 1 kink mode in a periodic cylindrical 
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Fig. 20. — Rotation frequency as a function of z 
in the nonlinear jet calculations discussed in Sec- 
tion [2] for various values of Vd, at times t = 
8.8,8.4,11.3 IV Rotation frequency is calculated 
by making a linear fit to the m = component of 
«6> for small values of r. 



plasma via linear MHD calculations. The linear 
calculations are treated as an initial value prob- 
lem in an Eulerian reference frame and as eigen- 
value problems in Eulerian and Lagrangian refer- 
ence frames. The results from all three methods 
are in agreement. While previous studies have 
shown that sheared flow is more efficient at sta- 



bilizing the kink mode (Wanex 2005), we show 



that rigid equilibrium rotation stabilizes the non- 
resonant m = 1 kink mode via the Coriolis effect. 
The Coriolis effect links radial and azimuthal mo- 
tions of the plasma, which distorts the kink eigen- 
mode and reduces its growth rate. 

The MHD equations used to model the jet prop- 
agation discussed in Section [2] include dissipative 
terms, and we should consider what effect dissipa- 
tion has on the rotational kink stabilization. In 
order to obtain smooth numerical solutions, the 
values chosen for the resistivity and the viscosity in 
the nonlinear jet simulations are much larger than 
that of any astrophysical jet system. However, 
we use the dissapationless ideal MHD equations 
for the eigenvalue analysis discussed in Section [4] 
While dissipation certainly affects the energy den- 
sities in the outflow in the jet simulations, the ro- 
tational stabilization is an ideal effect and robust 
to the choice of the dissipation coefficients. 

Our choice of initial conditions in the nonlinear 




Fig. 21. — The m = 1 Fourier component of v r in 
the nonlinear jet simulations with Vd = 0.5 and 
1.0 at times t — 8.76 and 10.26 T a respectively. 

jet simulations discussed in Sec. [2] has a signif- 
icant effect on the shape of the magnetic pitch 
profile, P(r), in the simulated jet. The combina- 
tion of inertia in the initial coronal plasma and a 
rapidly decreasing magnetic field acts as a back- 
ground which the magnetic flux can push against. 
This allows for the buildup of a large Bg , produc- 
ing a monotonically decreasing P(r). In contrast, 



the simulations of Moll ( 2008 ) produce jets with a 



monotonically increasing P(r). While these differ- 
ences affect the shape of the linear eigenfunctions, 
the eigenvalue calculations discussed in Sec. |4.2| 
show that the rotational stabilization is insensi- 
tive to the shape of the P(r) profile. 

With a decreasing P(r) profile and no equilib- 
rium rotation, there are lower and upper bounds 
on unstable values of k for the kink mode, and 
the growth rate, 7(A)), is a function of k. This can 
have a profound effect on the evolution of an ex- 
panding jet. The linear rigid rotation calculations 
discussed in this paper apply only to static equilib- 
ria. However, we contend that the results of these 
calculations can be used as a guide for considering 
the stability of the time-dependent equilibrium of 
an expanding jet. As the jet expands, the fc-value 
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of any given harmonic decreases in time, i.e. the 
harmonic is stretched by the jet expansion. If we 
consider an equilibrium that is expanding at a con- 
stant rate s with an initial length L and a mode 
with k = k' at time t = 0, the total energy gained 
by the harmonic over time t is can be estimated 
as 



AE(t,k') = E' I i 
la 



27(1^ 



(27) 



where E' is some initial energy in the mode. As 
we increase s, i.e. with faster expansion, the mode 
spends less time in the unstable range of k, and 
AE decreases. Moreover, equilibrium rotation 
acts to decrease the area under the "f{k) curve, 
decreasing AE as well. Clearly, this a nonlinear 
process, and the qualitative description given here 
motivates further study. 

While current driven instabilities may play a 
role in the wiggled structures which are observed 
in some outflows flReipurth |2002| |WorrarT"||2007| ); 
other explanations for these structures have been 
presented, such as precession of the the source ob- 
ject (Mascia dri ||2002 1 . In general, a combination 
of these effects could contribute to the formation of 
these structures. Since rotation is shown to stabi- 
lize the kink mode, knowledge of the jet rotation 
velocity relative to the Alfven velocity is critical 
for understanding the degree to which the kink 
plays a role. 
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